class: center, middle, inverse, title-slide # Minicurso ### Ana BVA, Rodolfo LCD ### 2021-03-22 --- <style type="text/css"> /* From https://github.com/yihui/xaringan/issues/147 */ .scroll-output { height: 80%; overflow-y: scroll; } /* https://stackoverflow.com/questions/50919104/horizontally-scrollable-output-on-xaringan-slides */ pre { max-width: 100%; overflow-x: scroll; } .remark-slide-content { font-size: 28px; padding: 20px 80px 20px 80px; } .remark-code, .remark-inline-code { background: #f0f0f0; } .remark-code { font-size: 24px; } .huge .remark-code { /*Change made here*/ font-size: 200% !important; } .tiny .remark-code { /*Change made here*/ font-size: 80% !important; } .center2 { margin: 0; position: absolute; top: 50%; left: 50%; -ms-transform: translate(-50%, -50%); transform: translate(-50%, -50%); } </style> --- # Aplicaciones <img src="./Images/RNA_seq_apps.png" width="400px" style="display: block; margin: auto;" /> --- # Expresión cuantitativa de loci (EQTL) -- - Relacionado con el estudio de rasgos y SNPs -- - Asociación entre cambios en la expresión de genes a SNPs conocidos -- - Determinar el efecto en el fenotipo -> mecanismos en *cis-* o en *trans-* -- - WGS y (variantes) y epigenómica (silencialmiento) para complementar los datos -- - Útil para la investigación de enfermedades o agronomía --  --- # Detección de variantes -- - Alternativa al estudio clásico de variantes (WGS o WES) -- - Presenta limitaciones relacionadas con la cobertura y la profundidad (cobertura mayor a 10X) -- - Se han reportado resultados en donde se identificaron > 8000 SNPs en músculo de *Longissimus thoraci* -- - Enfocado en la porción transcrita del genoma  --- # Redes de co-expresión -- - Asociación entre dos o más genes por medio de expresión coordinada -- - Genes co-expresados -> funciones relacionadas, mecanismos de regulación similar -- - RNA-seq permite obtener redes de co-expresión más completas que los microarreglos -- - Requiere amplio número de muestras (al menos 20) --- # Fusión de genes -- - Asociado a eventos de rearreglos genómicos como amplificaciones, deleciones y translocaciones -- - Región 5' gen A -> Region 3' gen B -- - RNA-seq permite conocer fusiones génicas en el desarrollo de enfermedades como cáncer -- - Permite conocer las diferencias en la expresión de genes fusionados entre dos condiciones  --- # Análisis de expresión diferencial - Aplicación más común del RNA-seq -- - Permite identificar y cuantificar la expresión de transcritos entre dos grupos distintos -- - Etapas del desarrollo, tejidos distintos, efecto de estímulos o tratamientos -- - Requiere del uso de modelos estadísticos que evaluen la expresión diferencial -- - Normalización de los datos (tamaño de la librería, longitud de los genes, composición de la librería...) --- class: inverse, center, middle # Tipos de archivos --- # Flujo de trabajo <img src="./Images/RNAseq_workflow.png" width="550px" style="display: block; margin: auto;" /> --- # FastQ -- - Derivado del formato FASTA -- - Id secuencia + **Calidad** -- - Generalmente no está comprimido (.fastq) > comprimir (.fastq.gz) -- - Comprende 4 líneas por secuencia: * @ ID del read + información de la corrida * Secuencia * Símbolo "+" * Calidad (Escala **Phred** y código **ASCII**) --- # FastQ ## Identificador -- - El identificador de la secuencia es la primera línea del archivo .fastq -- - Plataformas de Illumina superiores a 1.8 contienen la siguiente información: -- **@machneID:run number:flowCell ID:lane:tile:read:control number:index** ``` ## @SRR12038075.1 D00635:426:CD19YANXX:5:2301:1410:2130/1 ``` --- # FastQ ## Calidad La calidad de cada base secuenciada, en la escala **Phred**, se calcula de acuerdo la siguiente ecuación: .content-box-blue[ `$$Q = -10 * log10(p)$$` ] -- Q = Calidad -- p = Probabilidad de que el *base call* sea incorrecto -- Valores de **p** cercanos a cero -> Alta calidad -- ¿Cuál es la calidad si el valor p del basecall es de 0.01? -- `$$Q = -10 * log10(0.01) = -10 * -2 = 20$$` --- # Fastq ## Calidad <img src="./Images/Fastqc.png" width="700px" style="display: block; margin: auto;" /> --- # FastQ ## Calidad -- - Para evitar colocar los valores numéricos de la calidad en el archivo .fastq y saturarlo de información -> **[ASCII](https://www.ascii-code.com/)** -- - Phred + 33 -> Plataformas de Illumina a partir de 1.8 (símbolos especiales + alfanuméricos) -- - Phred + 64 -> Plataformas de Illumina anteriores a 1.8 <img src="./Images/ASCII0.png" width="400px" style="display: block; margin: auto;" /> <img src="./Images/ASCII33.png" width="400px" style="display: block; margin: auto;" /> --- # FastQ <img src="./Images/Fastq.png" width="700px" style="display: block; margin: auto;" /> --- # SAM/BAM -- - De manera posterior al alineamiento se generan archivos **"Sequence Alignment Map"** -> **SAM** -- - Archivos delimitados por tabuladores (TAB) -- - Convertidos a valores binarios como **"Binary Alignment Map"** -> **BAM** -- - **BAM** -> Mayor espacio en disco -> Acceso aleatorio más eficiente -- - Las secuencias en formato BAM se pueden ordenar por coordenadas -- - Conversion entre formatos [SamTools!](https://github.com/samtools/samtools) --- # SAM/BAM - Los archivos en formato **SAM** se dividen en dos secciones: -- - **HEADER**: Información general del archivo -- - **ALIGNMENT**: Información relacionada con el alineamiento de los reads -- <img src="./Images/SAM-BAM.png" width="550px" style="display: block; margin: auto;" /> --- <img src="./Images/Header.png" width="600px" style="display: block; margin: auto;" /> -- - La sección HEADER incuye unformación sobre cómo se generó el alineamiento y su almacenamiento -- - **Todas** las líneas del HEADER están separadas por tabuladores -- - La estructura general de las líneas es la siguiente: @ TAG:VALUE * TAG: Código de dos letras -> Contenido * VALUE: Valores del TAG --- # SAM/BAM ##HEADER Ejemplos de diferentes **TAG**: -- - @HD -> Metadatos del alineamiento -- - @SQ -> Datos de la secuencia de referencia -- - @RG -> Grupo de lecturas -- - @PG -> Información del alineador y el código del alineamiento --- <img src="./Images/Alignment.png" width="600px" style="display: block; margin: auto;" /> -- Esta sección procede al **HEADER** y contiene: -- - Cada línea corresponde a una lectura alineada -- - Para cada lectura se muestran los resultados de los siguientes campos obligatorios: <img src="./Images/Alignment2.png" width="400px" style="display: block; margin: auto;" /> --- # SAM/BAM ## FLAG -- - Los valores del FLAG contestan preguntas relacionadas con el alineamiento del read -- - Sí = 1, No = 0 -- - Cada respuesta se encuentra codificada en valores **binarios, decimales y hexadecimales** -- - Información útil para lecturas pareadas -> ¿La lectura es pareada? ¿Ambos reads alinearon adecuadamente? --- # SAM/BAM ## FLAG <img src="./Images/flag.png" width="600px" style="display: block; margin: auto;" /> --- # SAM/BAM ## Concise Idiosyncratic Gapped Alignment Report (CIGAR) - Almacena información relacionada con las operaciones necesarias para alinear el read - Se reporta con valores alfanuméricos: -- * M -> Alineamiento correcto * I -> Inserción comparado con la referencia * D -> Deleción comparado con la referencia * N -> Región con corrimiento (intrón) * S -> Soft clipping * H -> Hard clipping * P -> Padding o mutación silenciosa -- --- # SAM/BAM ## CIGAR <img src="./Images/cigar.png" width="600px" style="display: block; margin: auto;" /> --- class: inverse, center, middle # Organización de un proyecto bioinformático --- # Crisis de reproducibilidad en la ciencia La reproducibilidad en la ciencia implica que dos o más personas (o grupos de trabajo) realicen el mismo experimento bajo condiciones experimentales similares y obtengan resultados similares. Para ello se requiere los siguientes elementos: -- - Descripción de materiales y métodos -- - Datos crudos o semi-crudos -- - Código de cómputo utilizado -- - Versión del software utilizado -- - Información suplementaria para realizar el análisis de los datos --- # Crisis de reproducibilidad en la ciencia Actuamente, existe una **crisis de reproducibilidad** en la ciencia -> Perseguir explicaciones incorrectas -> Diseñar tratamientos incorrectos para ciertas enfermedades Para combatir este problema: -- - Transparencia de datos en publicaciones científicas -- - Apertura de la información (Open science) -- - Código abierto para el análisis de datos -- - Repositiorios o bases de datos públicas donde se almacene la información -- - Software depositado en repositorios o servidores a largo plazo --- # Crisis de la reproducibilidad en la ciencia <img src="./Images/Repo.png" width="400px" style="display: block; margin: auto;" /> Proyecto organizado en directorios y subdirectorios para reproducir el análisis de datos --- # Organización del proyecto ## README Un **README** es un archivo en texto plano (formado .md) que permite documentar el proyecto -- - El contenido de un proyecto (repositorio) -- - Las funciones de los scripts contenidos en el proyecto -- - El orden en el cual deben ser ejecutados los scripts para realizar los análisis -- - Software necesario para los análisis -- - Información suplementaria sobre los datos/scripts/software -- La función del **README** es explicar la organización del proyecto y que otros entiendan su contenido --- # Organización del proyecto ## README Tips: -- - Importante actualizar el README durante el desarrollo del proyecto -- - Editarlo en formato plano (Nano, TextWrangler, MacDown, ...) -- - Información clara y concisa <img src="./Images/README2.png" width="600px" style="display: block; margin: auto;" /> --- # Organización del proyecto ## Datos Los datos deben de almacenarse en un directorio especial (considerar espacio en disco): -- - Datos crudos -- - Datos semi-crudos -- - Datos procesados -- Considerar almacenar los datos en un repositorio como respaldo [OSF](https://osf.io/) --- # Organización del proyecto ## Datos Es recomendable llevar una tabla de metadatos durante el desarrollo del proyecto: -- - Almacenar información de cada muestra -- - Información de la corrida -> Lotes -- - Comentarios respecto a la corrida -- - Fechas de creación/modificación de los datos --- # Organización del proyecto ## Código **El código es una de las partes más importantes del proyecto** -- - Instrucciones que recibe el equipo para el análisis de datos -- - Diseñado a partir del pipeline o estrategia para abordar el análisis de datos -- - Elemento dinámico dentro del proyecto -- <img src="./Images/Script.png" width="400px" style="display: block; margin: auto;" /> --- # Organización del proyecto ## Código -- - Enumerar los scripts de acuerdo al orden en que serán ejecutados ``` +-- 1.QualityControl.sh +-- 2.Trimming.sh +-- 3.ReadAlignment.sh +-- 4.DifferentialExpression.R ``` -- - Utilizar rutas relativas hacia los archivos *input* ``` fastqc ../Data/*.fastqc -o ../Results/ ``` -- - Comentar el script -> facilita que sea entendido por humanos ``` #Change the location of bam files mv ../Data/*.bam ../Results/ ``` -- - ¡Mejorar el código! -> *Less is more* --- .center[] --- # Organización del proyecto ## Código .pull-left[ - **[DRY](https://en.wikipedia.org/wiki/Don%27t_repeat_yourself)** Don't Repeat Yourself - **[KISS](https://en.wikipedia.org/wiki/KISS_principle)** Keep It Simple Soldier - **[YAGNI](https://en.wikipedia.org/wiki/You_aren%27t_gonna_need_it)** You Aren't Gonna Need It - **[POLA](https://en.wikipedia.org/wiki/Principle_of_least_astonishment)** Principle of Least Astonishment - **[POLP](https://digitalguardian.com/blog/what-principle-least-privilege-polp-best-practice-information-security-and-compliance)** Principle of Least Privilege ] .pull-right[ <img src="./Images/Development.png" width="200px" style="display: block; margin: auto 0 auto auto;" /> ] --- # Organización del proyecto ## Resultados <img src="./Images/Results.png" width="700px" style="display: block; margin: auto;" /> --- # Organización del proyecto <img src="./Images/Proyecto.png" width="450px" style="display: block; margin: auto;" /> --- # Experimento Número de acceso GEO: GSE152699 --- class: inverse, center, middle # Alineamiento --- <img src="./Images/RNAseq_workflow2.png" width="600px" style="display: block; margin: auto;" /> Imagen tomada de [aquí](https://biocorecrg.github.io/RNAseq_course_2019/salmon.html) --- # Alineamiento <img src="./Images/ReadAlignment.jpg" width="600px" style="display: block; margin: auto;" /> Buscar la región del genoma/transcriptoma a partir de la cual se originaron los reads --- # Alineamiento - Indexar genoma/transcriptoma (FASTA, GTF3/GFF3) <img src="./Images/GenomeIndexing.png" width="500px" style="display: block; margin: auto;" /> -- - Alinear los reads (Fastq) al genoma/transcriptoma indexado <img src="./Images/ReadAlignment1.png" width="500px" style="display: block; margin: auto;" /> Imágenes tomadas de [aquí](https://almob.biomedcentral.com/articles/10.1186/s13015-016-0069-5/figures/1) y [aquí](https://figshare.com/articles/figure/_Comparison_of_stranded_and_unstranded_RNA_seq_library_methods_and_their_influence_on_interpretation_and_analysis_/1504417/1) --- # Alineamiento ## Spliced reads <img src="./Images/SplicedReads2.png" width="600px" style="display: block; margin: auto;" /> --- # Alineamiento Alineadores tipo splice aware -> Alineamiento de secuencias exon spanning -> Segmentan las secuencias en las regiones de dos exones - TopHat 2 <img src="./Images/SpliceAware.png" width="600px" style="display: block; margin: auto;" /> -- - HISAT <img src="./Images/HISAT.png" width="600px" style="display: block; margin: auto;" /> -- - STAR <img src="./Images/STAR.png" width="400px" style="display: block; margin: auto;" /> Imágenes tomadas de [aquí]() --- # Alineamiento ## Pseudo-alineadores Pseudo-alineadores (quasi-alineadores): - Kallisto - Sailfish - **Salmon** --- # Alineamiento ## Pseudo-alineadores <img src="./Images/Pseudoalignment.png" width="600px" style="display: block; margin: auto;" /> --- # Alineamiento ## Pseudo-alineadores <img src="./Images/PseudoAlignment1.png" width="500px" style="display: block; margin: auto;" /> -- <img src="./Images/PseudoAlignment2.png" width="500px" style="display: block; margin: auto;" /> --- # Alineamiento ## Pseudo-alineadores <img src="./Images/PseudoAlignment3.png" width="600px" style="display: block; margin: auto;" /> --- background-size: cover class: center, bottom, inverse background-image: <img src="./Images/hands_on.png" width="1315" /> --- # Hands on ## Generar índice del transcriptoma Para ello vamos a requerir: -- - Archivo **Fasta** del transcriptoma de referencia del humano. Descargado del sitio de [GeneCode](https://www.gencodegenes.org/human/) Lo pueden encontrar en el folder de `Transcriptome/` -- - Código para genear el índice ``` salmon index -t --genecode path_to_transcriptome.fa.gz -i path_to_save_index ``` `-t`: Ubicación al archivo Fasta del transcriptoma `-i`: Ubicación para salvar el índice `--genecode`: El Fasta del transcriptoma de referencia está en formato de GENECODE --- # Hands on ## Cuantificar la abundancia de los transcritos Requerimientos: -- - Archivos **Fastq** de las lecturas -> Paired end R1 y R2 Ubicados en el folder de `data` -- - Índice generado en el paso anterior -- - Código para producir cuantificar los transcritos ``` salmon quant -i path_to_index -l A -1 path_to_R1 -2 path_to_R2 -o path_to_store_results ``` `-i`: Ubicación del índice `-l`: Tipo de librería `-1 y -2`: Ubicación a las lecturas R1 y R2 `-o`: Ubicación para almacenar los resultados --- ```bash ##Crear un directorio para almacenar los datos de los conteos mkdir -p ../salmon_quant ##Llamar salmon para realizar el conteo salmon quant -i ../transcriptome/genecode.v37.salmon141 \ -l A \ -1 ../Data/s1_R1.fastq.gz -2 ../Data/s1_R2.fastq.gz \ -p 6 --validateMappings \ -o ../salmon_quant/s1_quant ``` --- class: inverse, center, middle # Objeto Summarized Experiment --- .center2[ <img src="./Images/Slots.png" width="700px" /> ] --- # Summarized Experiment <img src="./Images/Summarized_experiment.png" width="400px" style="display: block; margin: auto;" /> --- # Summarized Experiment <img src="./Images/ColData.png" width="400px" style="display: block; margin: auto;" /> --- # Summarized Experiment ## ColData -- - El cajón o slot correspondiente a `ColData` contiene la tabla de metadatos creada para importar las cuentas con tximeta -- - Para acceder a **.red[ColData]** usar el siguiente comando: .tiny[ ```r colData(se) ``` ``` ## DataFrame with 12 rows and 3 columns ## names condition cell ## <factor> <factor> <factor> ## A_Control_1 A_Control_1 Control A ## A_Control_2 A_Control_2 Control A ## A_Control_3 A_Control_3 Control A ## M_Control_1 M_Control_1 Control M ## M_Control_2 M_Control_2 Control M ## ... ... ... ... ## A_Verafinib_2 A_Verafinib_2 Verafinib A ## A_Verafinib_3 A_Verafinib_3 Verafinib A ## M_Verafinib_1 M_Verafinib_1 Verafinib M ## M_Verafinib_2 M_Verafinib_2 Verafinib M ## M_Verafinib_3 M_Verafinib_3 Verafinib M ``` ] --- # Summarized Experiment ## ColData - El slot **ColData** es un objeto con clase de *DataFrame* .tiny[ ```r class(colData(se)) ``` ``` ## [1] "DFrame" ## attr(,"package") ## [1] "S4Vectors" ``` ] -- - *Rownames* de **ColData** corresponden a los *Colnames* en el slot **Assay**- **.red[Importante para análisis con DESeq2]** .tiny[ ```r rownames(colData(se)) ``` ``` ## [1] "A_Control_1" "A_Control_2" "A_Control_3" "M_Control_1" ## [5] "M_Control_2" "M_Control_3" "A_Verafinib_1" "A_Verafinib_2" ## [9] "A_Verafinib_3" "M_Verafinib_1" "M_Verafinib_2" "M_Verafinib_3" ``` ] --- # Summarized experiment <img src="./Images/RowRanges.png" width="400px" style="display: block; margin: auto;" /> --- # Summarized experiment ## rowRanges - El cajón de `rowRanges` hace referencia a las coordenadas de cada transcrito -- - Para acceder al **.red[rowRanges]** usar: .tiny[ ```r rowRanges(se) ``` ``` ## GRanges object with 184844 ranges and 9 metadata columns: ## seqnames ranges strand | ## <Rle> <IRanges> <Rle> | ## ENST00000631435 CHR_HSCHR7_2_CTG6 142847306-142847317 + | ## ENST00000415118 14 22438547-22438554 + | ## ENST00000434970 14 22439007-22439015 + | ## ENST00000448914 14 22449113-22449125 + | ## ENST00000632524 CHR_HSCHR14_3_CTG1 105866322-105866332 - | ## ... ... ... ... . ## ENST00000444456 1 247974741-247975653 + | ## ENST00000416377 1 248003129-248004068 + | ## ENST00000318104 1 248121932-248122882 + | ## ENST00000606902 1 248498308-248498649 + | ## ENST00000427827 1 248549444-248549762 - | ## tx_id tx_biotype tx_cds_seq_start ## <character> <character> <integer> ## ENST00000631435 ENST00000631435 TR_D_gene 142847306 ## ENST00000415118 ENST00000415118 TR_D_gene 22438547 ## ENST00000434970 ENST00000434970 TR_D_gene 22439007 ## ENST00000448914 ENST00000448914 TR_D_gene 22449113 ## ENST00000632524 ENST00000632524 IG_D_gene 105866322 ## ... ... ... ... ## ENST00000444456 ENST00000444456 unprocessed_pseudogene <NA> ## ENST00000416377 ENST00000416377 unprocessed_pseudogene <NA> ## ENST00000318104 ENST00000318104 unprocessed_pseudogene <NA> ## ENST00000606902 ENST00000606902 unprocessed_pseudogene <NA> ## ENST00000427827 ENST00000427827 unprocessed_pseudogene <NA> ## tx_cds_seq_end gene_id tx_support_level ## <integer> <character> <integer> ## ENST00000631435 142847317 ENSG00000282253 <NA> ## ENST00000415118 22438554 ENSG00000223997 <NA> ## ENST00000434970 22439015 ENSG00000237235 <NA> ## ENST00000448914 22449125 ENSG00000228985 <NA> ## ENST00000632524 105866332 ENSG00000282455 <NA> ## ... ... ... ... ## ENST00000444456 <NA> ENSG00000237492 <NA> ## ENST00000416377 <NA> ENSG00000232215 <NA> ## ENST00000318104 <NA> ENSG00000177233 <NA> ## ENST00000606902 <NA> ENSG00000271934 <NA> ## ENST00000427827 <NA> ENSG00000227102 <NA> ## tx_id_version gc_content tx_name ## <character> <numeric> <character> ## ENST00000631435 ENST00000631435.1 83.3333 ENST00000631435 ## ENST00000415118 ENST00000415118.1 25.0000 ENST00000415118 ## ENST00000434970 ENST00000434970.2 55.5556 ENST00000434970 ## ENST00000448914 ENST00000448914.1 61.5385 ENST00000448914 ## ENST00000632524 ENST00000632524.1 54.5455 ENST00000632524 ## ... ... ... ... ## ENST00000444456 ENST00000444456.1 44.0307 ENST00000444456 ## ENST00000416377 ENST00000416377.2 45.3191 ENST00000416377 ## ENST00000318104 ENST00000318104.2 47.7392 ENST00000318104 ## ENST00000606902 ENST00000606902.2 45.9064 ENST00000606902 ## ENST00000427827 ENST00000427827.2 45.4545 ENST00000427827 ## ------- ## seqinfo: 455 sequences from GRCh38 genome ``` ] --- # Summarized experiment <img src="./Images/Assay.png" width="400px" style="display: block; margin: auto;" /> --- # Summarized experiment -- - El slot `assay` almacena la información de las cuentas para cada transcrito dividida en tres niveles: -- .tiny[ ```r assayNames(se) ``` ``` ## [1] "counts" "abundance" "length" ``` ] -- - Para acceder a la matriz de cuentas estimadas por *Salmon*, correr: .tiny[ ```r head(assay(se), 5) ``` ``` ## A_Control_1 A_Control_2 A_Control_3 M_Control_1 M_Control_2 ## ENST00000631435 0 0 0 0 0 ## ENST00000415118 0 0 0 0 0 ## ENST00000434970 0 0 0 0 0 ## ENST00000448914 0 0 0 0 0 ## ENST00000632524 0 0 0 0 0 ## M_Control_3 A_Verafinib_1 A_Verafinib_2 A_Verafinib_3 ## ENST00000631435 0 0 0 0 ## ENST00000415118 0 0 0 0 ## ENST00000434970 0 0 0 0 ## ENST00000448914 0 0 0 0 ## ENST00000632524 0 0 0 0 ## M_Verafinib_1 M_Verafinib_2 M_Verafinib_3 ## ENST00000631435 0 0 0 ## ENST00000415118 0 0 0 ## ENST00000434970 0 0 0 ## ENST00000448914 0 0 0 ## ENST00000632524 0 0 0 ``` ] --- # Summarized experiment ## Assay - Las matriz de abundancia *(TPM)* puede obtenerse: .tiny[ ```r ## Obtener matriz de TPM head(se@assays@data$abundance, 5) ``` ``` ## A_Control_1 A_Control_2 A_Control_3 M_Control_1 M_Control_2 ## ENST00000631435 0 0 0 0 0 ## ENST00000415118 0 0 0 0 0 ## ENST00000434970 0 0 0 0 0 ## ENST00000448914 0 0 0 0 0 ## ENST00000632524 0 0 0 0 0 ## M_Control_3 A_Verafinib_1 A_Verafinib_2 A_Verafinib_3 ## ENST00000631435 0 0 0 0 ## ENST00000415118 0 0 0 0 ## ENST00000434970 0 0 0 0 ## ENST00000448914 0 0 0 0 ## ENST00000632524 0 0 0 0 ## M_Verafinib_1 M_Verafinib_2 M_Verafinib_3 ## ENST00000631435 0 0 0 ## ENST00000415118 0 0 0 ## ENST00000434970 0 0 0 ## ENST00000448914 0 0 0 ## ENST00000632524 0 0 0 ``` ] --- # Summarized experiment ¿Recuerdan a qué tiene que ser igual *Rownames* del slot colData? -- .tiny[ ```r rownames(colData(se)) ``` ``` ## [1] "A_Control_1" "A_Control_2" "A_Control_3" "M_Control_1" ## [5] "M_Control_2" "M_Control_3" "A_Verafinib_1" "A_Verafinib_2" ## [9] "A_Verafinib_3" "M_Verafinib_1" "M_Verafinib_2" "M_Verafinib_3" ``` -- ```r colnames(assay(se)) ``` ``` ## [1] "A_Control_1" "A_Control_2" "A_Control_3" "M_Control_1" ## [5] "M_Control_2" "M_Control_3" "A_Verafinib_1" "A_Verafinib_2" ## [9] "A_Verafinib_3" "M_Verafinib_1" "M_Verafinib_2" "M_Verafinib_3" ``` -- ```r ## Comprobar que rownames de colData es igual a colnames de assay row.names(colData(se)) == colnames(assay(se)) ``` ``` ## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE ``` ] --- class: inverse, center, middle # Exploración de los datos --- .center2[ <img src="./Images/PCAplot.png" width="700px" /> ] --- # Exploración de datos ## ¿Por qué es importante explorar los datos? -- - Paso previo al análisis de expresión diferencial -- - Análisis de calidad de los datos -- - Permite conocer la congurencian entre individuos o réplicas -- - Mediante gráficos, visualizar comportamiento de los datos -> Outliers --- # Exploración de datos .pull-left[ <img src="./Images/PCAplot2.png" width="400px" /> ] -- .pull-right[ <img src="./Images/Heatmap.png" width="400px" /> ] --- # Exploración de datos .pull-left[ <img src="./Images/PCAplot2.png" width="400px" /> **.red[PCA]** ] .pull-right[ <img src="./Images/Heatmap.png" width="400px" /> **.red[Heatmap]** ] --- # PCA -- - Análisis de componentes principales -- - Método algebraico para reducir la dimensionalidad de sets de datos complejos (múltiples variables) -- - Reducción de dimensionalidad o variables permite un análisis exploratorio más intuitivo -- - Reducción de variables implica preservar y captar la mayor información sobre los datos --- # PCA ## Normalización de los datos -- <img src="./Images/zscore.png" width="350px" style="display: block; margin: auto;" /> Escalación de los datos -> Normalización para hacerlos más comparables -- <img src="./Images/covariancematrix.png" width="350px" style="display: block; margin: auto;" /> Busqueda de correlación entre las variables -> correlación positiva o negativa --- # PCA ## Proyección de los datos .center[  ] -- .content-box-red[Los componentes principales son los nuevos ejes que maximizan la distancia de los datos al origen ] --- # PCA ## Componentes principales -- - El .green[primer componente] es aquel en el que los datos presentan la mayor separación (variación) -- - El .orange[segundo componente] es el que tiene la segunda mayor separación entre los datos y es perpendicular al primer componente -- - ¿Cuantos componentes existen? -> tantas variables en el set de datos --- # PCA ## Componentes principales - Cada componente principal tiene asociado un eigenvector (vector unitario) y un eigenvalue (cantidad escalar) -- - Los eigenvalues son la suma del cuadrado de las distancias de los puntos proyectados sobre dicho componente principal -- - Los eigenvalues permiten seleccionar cuál es el componente principal que explica la mayor variación en los datos --- # PCA <img src="./Images/Screeplot.png" width="450px" style="display: block; margin: auto;" /> Screeplot --- # PCA - En sets de datos de RNAseq, el PCA permite visualizar la distancia o congruencia de los datos -- - Generalmente son gráficas de puntos (muestras) en dos dimensiones (dos componentes) que resumen las principales fuentes de varianza -- .center[ <img src="./Images/PCARNAseq.png" width="400px" /> ] --- # PCA ## ¿Cómo crear nuestro propio PCA? - Primero instalen las siguientes librerías: .tiny[ ```r library(DESeq2) library(PCAtools) ``` ] -- - Usemos la siguiente normalización logarítmica de las cuentas (toma en cuenta el tamaño de la librería) .tiny[ ```r ## Transformación recomendada para sets de con menos de 30 muestras rld <- rlog(dds, blind = F) ``` ] - Si tuvieramos un número mayor de muestras (n > 30) entonces conviene: .tiny[ ```r vsd <- vst(dds, blind = T) ``` ] --- # PCA - Usemos la función interna de *DESeq2* para graficar nuestro PCA: .tiny[ ```r ## El argumento intgroup permite especificar mediante cuál variable colorear los datos plotPCA(rld, intgroup = "condition") ``` ] .content-box-red[*DESeq2* por medio de `plotPCA` evalua la varianza de los 500 genes con mayor varianza del set de datos ] - Con la librería de *PCAtools*: .tiny[ ```r ## Crear un objeto que contenga los datos del PCA PCA <- pca(assay(rld), scale = T, metadata = coldata) ## Graficar en 2D los resultados biplot(PCA, colby = "condition") ## Generar un screeplot para visualizar la varianza asociada a cada componente screeplot(PCA) ``` ] --- .center[ ] --- # PCA ## ¿Cómo ven los resultados? -- .pull-left[**DESeq2** <img src="Slides_files/figure-html/unnamed-chunk-66-1.png" width="350px" style="display: block; margin: auto;" /> ] -- .pull-right[**PCA tools** <img src="Slides_files/figure-html/unnamed-chunk-67-1.png" width="350px" style="display: block; margin: auto;" /> ] --- # PCA ## ¿Cuánta varianza es explicada por cada componente? <img src="Slides_files/figure-html/unnamed-chunk-68-1.png" width="400px" style="display: block; margin: auto;" /> --- # PCA y MDS .scroll-output[ Generemos un MDS plot interactivo con los datos. Para ello usemos la librería de *Glimma* y la función de `GlimmaMDS` ```r glimmaMDS(dds) ```
] --- class: inverse, center, middle # Visualización de los resultados --- .center[ ### ¿Cuántos tipos de gráficos conoces para visualizar resultados del análisis DE? ] .pull-left[ <img src="./Images/MAplot.png" width="350px" /> ] .pull-right[ <img src="./Images/Volcanoplot.png" width="350px" /> ] --- .center[ ### ¿Cuántos tipos de gráficos conoces para visualizar resultados del análisis DE? ] .pull-left[ <img src="./Images/MAplot.png" width="350px" /> **.green[MAplot]** ] .pull-right[ <img src="./Images/Volcanoplot.png" width="350px" /> **.green[Volcano plot]** ] --- # MAplot - Gráfico que representa la distribución de los genes o transcritos en las comparaciones de interés - **M** eje y de *minus* -> `$$logTx - logCT = logTx/CT$$` - **A** eje x de *average* -> Promedio de las cuentas normalizadas para cada gen en todas las muestras - Para generar el **MAplot** usemos el siguiente comando: .tiny[ ```r ##Es importante que recuerden que la hipótesis nula que se probó fue ##"El lfc del gen n es igual a 0" por lo tanto los genes coloreados son... plotMA(res) ``` ] --- # MAplot <img src="Slides_files/figure-html/unnamed-chunk-75-1.png" width="450px" style="display: block; margin: auto;" /> --- # MAplot Para tener una mejor estimación del LFC de genes que: - Son poco abundantes - Tienen alta dispersión Usemos la función `lfcShrink` de la librería *.red[apeglm]* .tiny[ ```r ##En el argumento coef indicar el contraste entre los grupos experimentales res <- lfcShrink(dds, coef = "condition_Verafinib_vs_Control", type = "apeglm") ##Si desconocemos el nombre del contraste, usar: resultsNames(dds) ``` ] --- # MAplot .pull-left[ <img src="Slides_files/figure-html/unnamed-chunk-77-1.png" width="450px" style="display: block; margin: auto;" /> ] .pull-right[ ```r plotMA(res2) ``` <img src="Slides_files/figure-html/unnamed-chunk-78-1.png" width="450px" style="display: block; margin: auto;" /> ] --- # MAplot .scroll-output[ También pueden generar un MAplot interactivo con la librería de *Glimma* ```r glimmaMA(dds) ``` ``` ## 343 of 36927 genes were filtered out in DESeq2 tests ```
] --- # Volcano plot De manera similar al MAplot con el volcano plot visualizamos los genes que muestran expresión diferencial en nuestra condición de interés - En el eje y se grafica el -log10 de padj - En el eje x se grafica el lfc o *log2foldchange* .tiny[ ```r ##Para crear un volcano plot necesitas convertir los resultados de DESeq a un data frame DEG <- as.data.frame(res) ##En el script de funciones encontrarás la función de volcanoplotR para generar tu gráfico #los valores de los argumentos logfc y p.adj deben ser iguales al threshold utilizado para generar los resultados #en type debes indicar que los resultados provienen de DESeq volcanoplotR(DEG, logfc = 0, p.adj = 0.1, type = "DESeq") ``` ] --- # Volcano plot <img src="Slides_files/figure-html/unnamed-chunk-81-1.png" width="450px" style="display: block; margin: auto;" /> --- # Heatmap El *heatmap* nos permite visualizar la expresión de los genes diferencialmente expresados en terminos de las cuentas normalizadas Consideraciones: - Usar los valores de las cuentas normalizadas para una mejor comparación entre muestras - Escalar los valores de las cuentas (renglones) para visualizar las diferencias en la expresión Usaremos la librería de *.red[pheatmap]* --- # Heatmap .tiny[ ```r ##Guardar la lista de transcritos que mostraron expresión diferencial significativa significant <- res_df %>% filter(log2FoldChange > 0 & padj < 0.1 | log2FoldChange < 0 & padj < 0.1) ##Cortar la matriz de cuentas normalizadas con los id de los transcritos diferencialmente expresados norm_counts <- norm_counts[rownames(significant)] ##Generar una tabla de anotaciones preservando el tratamiento y el tipo de células annotation_col <- coldata[, c(3, 4)] ##Generar el heatmap empleando clustering jerarquico pheatmap(norm_counts, border_color = NA, scale = "row", clustering_distance_rows = "euclidean", clustering_distance_cols = "euclidean", clustering_method = "average", show_colnames = F, show_rownames = F, annotation_col = annotation_col) ``` ] --- # Heatmap <img src="Slides_files/figure-html/unnamed-chunk-83-1.png" width="450px" style="display: block; margin: auto;" />